Exact field ionization rates in the barrier suppression-regime 
from numerical TDSE calculations 



Numerically determined ionization rates for the field ionization of atomic 
hydrogen in strong and short laser pulses are presented. The laser pulse in- 
tensity reaches the so-called "barrier suppression ionization" regime where 
field ionization occurs within a few half laser cycles. Comparison of our nu- 
merical results with analytical theories frequently used shows poor agree- 
ment. An empirical formula for the "barrier suppression ionization" -rate is 
presented. This rate reproduces very well the course of the numerically deter- 
mined ground state populations for laser pulses with different length, shape, 
amplitude, and frequency. 
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I. INTRODUCTION 



With the "table-top" laser systems, nowadays available, laser pulse peak field strengths 
much greater than the binding field of the outer atomic electrons can be achieved (see e.g. |l ] 
for an overview). Above a certain threshold electric field the electron is able to escape even 
classically from the atomic nucleus, i.e., without tunneling through the barrier formed by 
the Coulomb potential and the external electric (laser) field. This regime is called "barrier 
suppression ionization" (BSI) [0. 

In combination with the dramatic progress in decreasing the pulse duration below 10 fs 
|3|-|| new features in the ionization dynamics are expected. In particular, ionization at such 
high field strengths occurs mainly within a few half laser cycles, i.e., on a sub-femtosecond 
time scale, provided that the pulse rises fast enough so that tunneling contributes negligibly 
to the overall ionization. Fast depletion of bound states within one half laser cycle leads to a 
non-isotropic electron distribution. Apart from the peaked angular distribution of the photo 
electrons in electric field direction, in the BSI-case there is also an asymmetry along this 
field axis 0. This opens up the possibility to manipulate the electron distribution function 
of laser produced plasmas. By "tailoring" the pulse shape the plasma formation process may 
be controlled according to the application under consideration, e.g., harmonics generation 
||, or XUV laser schemes ||. 

Experimentally observed ion yields are usually analyzed by means of tunneling theo- 
ries among these Ammosov-Delone-Krainov (ADK) flQfl , Keldysh ||11|| , Keldysh- Faisal- Reiss 
(KFR) jy! or Landau [13|] theory are the most prominent ones. However, it is, in general, 
not possible to get good agreement for several ion species without "shifting" the laser in- 
tensity 0. By examining the derivations of KFR-type theories it becomes obvious that 
they should fail in the "barrier suppression ionization" (BSI) regime because the transition 
between an unperturbed initial state and a Volkov state is calculated there. However, the 
influence of the strong laser field on the inneratomic dynamics must not be neglected in BSI. 



An attempt to extend the ADK-theory to BSI has been undertaken A pure classical 
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ionization rate has been proposed recently |T5 . 

In this paper we compare numerically determined ionization rates for various kinds of 
pulse shapes and peak field strengths with results predicted by several analytical derivations: 
the Landau tunneling formula [[HJ, the Keldysh rate |ITJ| , the ADK formula JIIJ and its 



extension to the BSI-regime JR], a classical rate derived by Posthumus et al. [T5|] and a 



tunneling rate suggested by Mulser fl6fl . In our numerical studies we restrict ourselves to 
the ionization of atomic hydrogen in an intense, short, linearly polarized laser pulse. We 
focus on the field strength region where the ionization rate is of the order of the laser 
frequency because ionization occurs within a few half laser cycles in this case. 

In Section || we review the time-dependent Schrodinger equation (TDSE) of field ion- 
ization. Moreover, we state the analytical formulas used for comparison with our numerical 



results. In Section |TTTj we present our numerical results for various pulse shapes and field 
strengths. The numerical results are discussed in Section [TV|. We conclude in Section [V]. 
Details on the numerical method are attached in the Appendix. 

II. THEORY 

A. Time-dependent Schrodinger equation (TDSE) 

The TDSE for an electron interacting with the nuclear potential —Z/r and the laser field 
E{t) in dipole approximation and length gauge reads [[L7j 

i^*(r,*)= f-^-| + rJ0(f)W,f) (1) 



(atomic units (a.u.) are used throughout this paper |]18| ). If the electric field is chosen 
to be directed along the z-axis, cylindrical coordinates are introduced and the Ansatz 
ty(p,(p,z,t) = ip(p, z, t) exp(im<^)(27r)~ 1//2 is made, the TDSE assumes the following two- 
dimensional form, 

■ 9 , 1 (I d , d n m 2 d 2 \ , . _ , Z . ( 



and the normalization condition 

/*oo 

dp p dz \ip{p, z, t)\ 2 = 1 (3) 



holds. The TDSE (|2|) was numerically solved first by Kulander in 1987, but for intensities 
below 10 15 W/cm 2 0. 



In a recent work by Kono et al. [^(J it was systematically examined for what parameter 
A the substitution 

^,z,t) = v^e~ i/2 ^\z,t), e=p, z= z , t=t, (4) 

is most favorable numerically. It turned out that the choice A = 3/2 is best, both for stability 
and accuracy. The TDSE corresponding to the substitution (^) is given in Appendix [A[ We 
used a Peaceman-Rachford scheme to propagate the wavefunction z,t) (see Appendix 
[A|or Ref. [OT for details). Absorbing boundary conditions were implemented which keep the 
main interaction region in the vicinity of the atomic nucleus free from otherwise reflected 
probability density. 

In all our calculations we started from the Is ground state, i.e., m = 0. The stable 
ground state on the numerical grid (which is slightly different from the analytical solution 
of the Coulomb problem, depending on the grid-spacing ) was determined by applying our 
propagation scheme with an imaginary timestep to the grid representation of the known 
analytic solution. 

B. Ionization rate formulas 

In this Section we review the ionization rate formulas used for comparison with our 



numerical results of Section 111. If we assume that an ionization rate W^-E^i)] is given, the 



probability for the electron to remain bound is 

r(t) = exp (- jT W[E{t')} dA . (5) 



We take 



A(t) = i - r(t) (6) 

as the ionization probability which is, apart from a small time-shift, equivalent to the com- 
mon procedure to calculate the amount of probability to find the electron in a small volume 
around the atomic nucleus, 

ra ra 

A'(t) = l- dp pi dz\^(p,z,t)\ 2 , a^5a.u. (7) 

JO J -a 

We assume that the laser pulse "hits" the atom at t = (or ionization is negligible for 
t < 0). 



1. Landau formula 

Landau & Lifshitz derived a formula for the ionization rate of hydrogen when the electron 
is in the ground state initially ||13| . The result is easily extended to hydrogen-like ions (where 
the ground state energy is £ = —Z 2 /2), 



2. Keldysh formula 

Keldysh perturbatively calculated the transition rate from an initial bound state to a 
state representing a free electron in a laser field ( Volkov state) |~T| , 

w (67r)1/ V ( E V /2 ( 2 ^\f l2 \ ^ 



3. Ammosov-Delone-Krainov (ADK) formula 

Ammosov, Delone, and Krainov derived a tunneling ionization rate for complex atoms 
in an ac electric field |10|| . The initial state is described by an effective quantum number 
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n* and the angular and magnetic quantum numbers I and m, respectively. The ADK-result 
reads 

W ADK = Cl. t W, mm (^pg^)" 2 (|(2|4l)^) 2 """ W " 1 exp (10) 
with 

'2e\ n 



/ ze \ 

Cn^=(-J (27m*)- 1 / 2 , /(/,m) 



(2£ + !)(£ + 


m\)\ 


0\m\ 


m 




m\)\ 



The constant e in the coefficient C n *i is Euler's number 2.71828... In the derivation of the 
ADK rate ([31]) averaging over one laser cycle was performed. The validity of the ADK- 
formula is expected to be best for n* 1, E -C 1, and uj <C \S \. 

4- BSI extension to ADK 
Krainov suggested an extension of ADK-theory to incorporate BSI [14]. The result is 

WK! = f^i £ °D 3/2 V"' r ^ U + JJMj) in, 

irn* (2E) 1 /3 ^ En* J Jo \ (2£) 3 / 2 J v ; 

where Ai denotes the Airy-function. Formula ( |TT]) reduces to the usual ADK rate (pi]) in 
the limit of a relatively weak laser field (tunneling limit). 

5. Classical rate proposed by Posthumus et al. 



Recently, Posthumus and co-workers proposed a purely classical BSI ionization rate []15 
Taking the equipotential surface corresponding to the atomic ground state and examining 
its intersection with the field-deformed Coulomb potential enables the authors to calculate 
the rate from a geometrical viewpoint. Their result reads 

_ 1 - £ 2 J(AZE) _ nZ 

Wd ~ 27L ' T ° " \£o\(2\£ \)^- (12) 

T is the classical orbit period for the so-called "free falling" -trajectories with zero angular 
momentum. The authors of [|15j present also a cycle-averaged expression of the rate. They 
finally suggest to take W c \ + W / adk(-^ci) as the total ionization rate in the BSI regime where 
J c i is an appropriate threshold intensity. 



6. Tunneling rate proposed by Mulser 



Mulser calculated the ionization rate by approximating the tunneling barrier formed 



by the Coulomb potential and the external field with a barrier parabolic in shape [16 
After calculating the transmission coefficient through this parabolic barrier and making an 
assumption for the tunneling current the rate formula 

\E 1 . A + exp|/3| ( 7 -3a \ . . 

^Mu = ^hi A+ \ , where A = exp ^ —C) , (13) 

3 + a 4E 1 / 2 > , cn i/g 2|goj 

P=— C > a= (2^7i' C=-^{2\£ \) ¥J ^ Ji , 

is obtained. 



III. NUMERICAL RESULTS 

In this Section we study the ionization dynamics of the Is atomic hydrogen electron 
under the influence of the external laser field E(t). The laser field is assumed to have the 
form 

E(t) = E(t) sin(cjt + <f) (14) 

where E(t) is the pulse shape-function and <p is a constant phase. In the following we vary 
the pulse envelope E(t), the laser frequency u, and the phase ip in order to examine their 
influence on the temporal evolution of the ground state probability T(t). 



A. Instantaneously switched on dc field 

Although the dc field instantaneously switched on is, from the experimental point of 
view, not realistic at all, this case delivers useful insight in how important transient effects 
might be. Furthermore it is interesting to check whether ionization occurs with a constant 
rate after transient effects have died out. 
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In the instantaneously switched on field-case the envelope function is 

E(t) = E = const. for t > (0 otherwise). (15) 

In Figure [I] the ground state population T(t) is plotted vs time for the five different 
amplitudes E = 0.1, 0.2, 0.3, 0.4 and 0.5. One easily verifies that after a very short transient 
period of about 2 a.u.= 0.048 fs the constant rate-behavior sets in. This transient time 
period may be estimated by purely classical considerations if one assumes that the atomic 
response time is similar to that of a classical system with an electron density corresponding 
to the quantum mechanical probability density of the ground state. The electron density 
then is n e rs (47r/3) _1 a.u. which leads to a "plasma frequency" cu p fa 3 1 / 2 a.u. The classical 
response time therefore would be about 3.6 a.u. = 0.09 fs. 

The constant rates W are given in the plot. We postpone a comparison with the analyt- 
ical rate formulas mentioned above till Section |TV|. 

The probability density \ip(£, z)\ 2 after 15.5 atomic time units for the E = 0.3-case is 
shown in Figure 0. Since we chose E > the electron escapes in negative z-direction. Note 
the pronounced asymmetry and the Is peak which does not move as a whole; it rather 
persists at the Coulomb singularity. 

B. Square pulses and phase dependence 

Now we study an ac field with a step-like envelope function, 

E(t) =Esm(ut + (p) for t>0 (0 otherwise). (16) 

In Figure |3] the ground state populations for the two field amplitudes E = 0.3 and 0.5 are 
shown. In each case three different phases (<p = 0, 7r/4, tt/2) were chosen in order to check 
how strong ionization depends on phase effects. The frequency was u = 0.2 in these runs. 

During the course of one half cycle ionization is strongly phase dependent. In the E(t) = 
E cos i^t-cases ionization is particularly strong in the beginning owing to the abrupt turn- 
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on of the field, while in the E(t) = E sin cut-cases ionization starts smoothly. A steady 
state-rate, based on cycle-averaging, of course cannot resolve such details. 

For E = 0.3 ionization lasts mainly two half cycles while for E = 0.5 already after one 
single half cycle ionization is > 98%. The more rapid ionization is, the stronger should be 
the dependence of ionization on the phase (p. However, even in the E = 0.5-case the two 
fields E(t) = sin cut and E(t) = -E 1 cos cut lead to the same net ionization after one half 
cycle. Only if one is interested in the ionization dynamics on time scales below one optical 
half cycle ionization becomes phase-sensitive. However, even the shortest pulses nowadays 
available have to cross the field strength region where the ionization rate is ~ cu. Once this 
regime is passed there is not much electron density left to be ionized by the stronger part 
of the pulse. 

For the sake of illustration the probability density after one complete optical cycle in the 
E(t) = 0.3 cos cut-case is shown in Figure |j. Owing to rescattering of probability density at 
the ionic core wave-packets have already built up. Closer examination yields that subsequent 
wave packets in position space can be mapped to subsequent wave packets in momentum 
space. These momentum space-packets differ in energy by the amount of hu, and thus are 
the famous "above threshold ionization" (ATI)-peaks (see |21[] for a detailed analysis). 



C. Gaussian pulses 

A shape which resembles in a reasonable manner an experimental laser pulse is Gaussian. 
We took 

E(t) = E(t) sin cut, E(t) = Eexp (- ^^ J . (17) 

Since a Gaussian is infinitely extended we have to start our computer runs with non- vanishing 
E(0). We chose E(0) to be 5% of the maximum field amplitude E. Demanding the Gaussian 
envelope to cover N laser cycles within the region E(t) > 0.05.E yields 

t = A%/cu, a 2 = tg/(41n20). (18) 



In Figure ^| the ground state populations for the four Gaussian pulses with E = 0.3, 0.5 
and N = 6, 12 each and to = 0.2 are shown. Besides, the result for a lower frequency 
(u = 0.1) and E = 0.5, N = 12 is included. The 12-cycle E = 0.3-pulse (drawn solid) 
ionizes most slowly, but the 6-cycle E = 0.3-pulse (dotted) deplets quicker the ground state 
than it is the case for the 12 cycle E = 0.5-case (dashed). This is due to the fact that the 
BSI regime is reached earlier for the weaker but shorter E = 0.3-pulse. 

The low frequency pulse (thin solid line) causes more rapid ionization than its counterpart 
with twice the frequency since the total time where the BSI region is reached (measured in 
absolute time units) is larger. 

We will further discuss the ground state populations depicted in Figure [| in Section [TV 
when we reproduce them with an empirical formula. 

D. Sin 2 pulses 

In Ref. one of the authors (D.B.) dealt extensively with sin 2 -pulses of the form 

E(t) = £sin 2 (^tj smut, T = Nx^. (19) 

Since the results look very similar to those in the Gaussian case we suppress a further 
discussion here. However, in Section [TV] we utilize rates numerically determined in J!]] for 
sin 2 -pulses in order to confirm the insensitivity of our proposed rate formula with respect 
to the pulse shape. Furthermore, a different numerical scheme was used in J7J. This gives 
additional reliability to the numerical results which will be utilized to derive an empirical 
BSI rate in the following Section. 

IV. DISCUSSION 

In this Section we want to demonstrate that it is possible to reproduce our numerical 
results using a simple formula for the ionization rate in the BSI regime. This rate is not 
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sensitive to laser frequency and pulse shape in a wide parameter range. Moreover we show 
that none of the analytical rates stated in Subsection [11 B| is applicable to BSI. 



Since BSI occurs mainly during one or two half laser periods a cycle-averaged rate ob- 
viously makes no sense. Therefore the laser field E(t) with its entire time- dependence has 
to be plugged in a rate formula, i.e., W(t) = W[E(t)}, while in tunneling ionization a rate 
which depends on the pulse envelope only, W{t) = W[E(t)}, is sufficient. 

We determined instantaneous ionization rates from the decreasing ground state popula- 
tions, in accordance with Eq. (|6]). In Figure |^ the results are plotted vs the electric field 
present at the corresponding instant. Usually the deepest descent in the ground state pop- 
ulation is in the vicinity of the electric field maximum of the actual half cycle. However, 
this behavior might be disturbed by "backsweeping" probability density ionized earlier, es- 
pecially for high frequencies (frequencies not much less than |£ |) since the excursion length 
of a freely oscillating electron is then not much larger than the width of its wave-packet 
representation. 

In Figure |^ different symbols are used for different pulse shapes, pulse lengths, and laser 
frequencies. For comparison the predictions by the analytic formulas of Subsection |11B 
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drawn as well. The scattering of the numerical data is due to the fact that instantaneous 
rates for a certain electric field value may stem from runs with different pulse shapes, peak 
field strengths or laser frequencies. 

The BSI regime for atomic hydrogen sets in for E = 0.146 when a classical electron, 
initially on an £ = —0.5 orbit, can escape from the atomic core. In general this so-called 
critical field in the case of hydrogen-like ions is given by |22|-p4 



E crit = (V2 + l)\£ \ 3 / 2 . (20) 

Once the critical field is reached one expects rapid ionization within a few half cycles. 
Therefore we are especially interested in the region where E > 0.15. Fortunately, the 
scattering of our numerical data is small in this region of field strengths. This makes possible 
our goal to provide a BSI rate formula valid for a wide range of pulse shapes and laser 
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frequencies. 



We observe that none of the analytical theories under consideration predicts the BSI rates 
correctly. Apart from Keldysh's result all formulas overestimate the ionization rate in the 
region of interest, 0.15 < E < 0.5. The ionization rate for much higher field strengths might 
be of academic interest since such high field strengths cannot be reached without strongly 
ionizing the hydrogen atom during earlier parts of the pulse where the field strength is in the 
region we focus on in this paper. In real experiments, with rare gases for instance, there are 
of course stronger bound electrons which get free not before E ^> 0.1 but for those electrons 
.Ecrit is larger too. We will discuss the scaling behavior of the ionization rate with respect 
to Z lateron. 

The rates of Posthumus (P) and Mulser (M) saturate at higher field strenghts. This 
is owing to taking the unperturbed inneratomic motion to derive an ionization current. In 
reality, however, the external field influences the inneratomic motion of the electron and 
yields a higher ionization current. The tunneling theories (L, Al and A2) are even worse 
when extrapolated to higher field strengths; they predict a decreasing ionization rate which 
is clearly unphysical. Note that "stabilization" cannot occur when ionization lasts less than 
one laser cycle. Although the Keldysh rate (K) does not suffer from these shortcomings it 
underestimates the ionization rate by a factor three and more. 

The numerically determined ionization rates in the region 0.15 < E < 0.5 can be nicely 
fitted by W = 2AE 2 . Since every realistic pulse passes through a region where the electric 
field is within the tunneling regime we propose a combined formula 



where E' is a threshold electric field determined by imposing W(t) to be continuous, and 
W'(t) is an appropriate tunneling rate. For the Landau rate E' = 0.084 holds. 



Landau tunneling rate. The agreement with the exact numerical results (drawn dotted) is 




(21) 



In Figure [7| the solid curves were calculated by applying the BSI rate ( |21"|) to the four 
Gaussian pulses which led to the results already depicted in Figure |5|. For W we used the 
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satisfactory. Deviations, especially in the N = 12, E = 0.3-run, are mainly due to the (even 
for lower field strength) not very accurate Landau rate. For shorter pulses and higher peak 
field strengths the agreement becomes excellent. The dashed curve is the result when the 
Landau rate alone is applied to the entire N = 6, E = 0.5 pulse; the ionization rate is 
strongly overestimated. 



In Figure || the BSI rate (|2lD was evaluated for the square pulses discussed in Subsection 



111 B[ In the upper plot the agreement with the numerical results for the E sin cut-case is 



good. However, in the lower plot [E cos cut-case) the agreement is not particularly good since 
the abrupt jump in the field strength from (for t < 0) to E for (t > 0) leads to transient 
dynamics which cannot be reproduced by our simple rate (^l|). Therefore, care has to be 
excercised for laser pulses where the BSI regime is reached rather abruptly on time scales 
shorter than one quarter laser cycle. In all other cases the rate formula (|2l|) worked well. 

A. Scaling 

The TDSE ([[]) can be rescaled to the atomic hydrogen-case by substituting 

r = Zr, t = Z 2 t, Co = cu/Z 2 , E = E/Z 3 . (22) 

Since our BSI rate is not sensitive to u, and an ionization rate has the dimension of an 
inverse time the rescaled result reads 

, W'[E(t)] for E(t) < E' 
W(t) = { L WJ W (23) 

2A/Z 4 x E{t) 2 for E(t) > E' 

B. The role of the Keldysh parameter 

The Keldysh parameter 

T -(Sf (24) 
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with Up the "ponderomotive potential" U v = E 2 /(4:Ui 2 ), i.e., the mean quiver energy of an 
electron in the laser field, has to be much less than unity when tunneling theories such as 
ADK are derived. The Keldysh parameter has the vivid physical interpretation of tunneling 
time measured in units of the laser period. Does the Keldysh parameter reveals some 
significance in the BSI regime too? First of all we note that the Keldysh parameter in our 
numerical examples is not much less than unity. In the E = 0.3, u = 0.2-case it is 0.67, 
in the E = 0.5, u = 0.1-case it is 0.2. Thus, in commonly used terms in this field, we are 
rather in the multiphoton than in the (to BSI extended) tunneling regime. 

However, the static field-rates in Fig. p] are also well covered by our empirical rate. 
Additional test runs at intermediate frequencies yielded good agreement also. Thus, the 
insensitivity of our BSI rate with respect to the laser frequency (from static fields up to 
uj = 0.2) shows that there seems to be no need to put much emphasis on the concept of the 
Keldysh parameter in BSI. However, we did not deal with frequencies > |£o| i n this paper. 
Moreover, a small laser frequency keeps the portion of already ionized probability density far 
away from the ionic nucleus most of the time since the excursion length is large. Therefore, 
the ionization curves for smaller frequencies usually look "cleaner" since interference with 
parts of the wave function representing the already ionized electron is suppressed. 

V. CONCLUSION 

We conclude that even for the simplest atom we can think of, i.e., atomic hydrogen, none 
of the theories discussed in this paper predict correctly the ionization rate in short intense 
laser pulses reaching the BSI regime. Thus, extrapolation of tunneling theories to BSI is not 
permitted. From the numerical results we deduce that a successful theory should take the 
influence of the strong laser field on the inneratomic dynamics into account. For quantum 
treatments of strong field ionization this means that one must not make the assumption that 
the initial state (to be plugged into the transition matrix element) evolves in time as if it 
was unperturbed (as it is usually done in KFR-type theories). In classical theories (such like 
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the one by Posthumus et al.) this corresponds to taking the effect of the laser field on the 
bound Kepler-orbits into acount. However, since in either case, quantum or classical, this 
appears extremely hard to achieve, empirical rates from numerical simulations of strong field 
ionization are highly desirable and important as an incredient for other simulation codes, 



In this paper, an empirical formula for the BSI rate has been proposed. Our formula 
is not sensitive to pulse shapes and laser frequencies in a wide parameter range, especially 
when combined with a reliable tunneling formula for the weaker parts of the laser pulse. 



This work was supported in part by the European Commission through the TMR Net- 
work SILASI (Super Intense Laser Pulse-Solid Interaction), No. ERBFMRX-CT96-0043 and 
by the Deutsche Forschungsgemeinschaft (DFG) under contract no. MU 682/3-1. 

APPENDIX A: NUMERICAL METHOD 

Starting point is the TDSE (|2|). We follow the line of Kono et al. ]2IJ and perform the 
substitution (f|) 
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Z = Z, t = t. 



(Al) 



The normalization condition for z, t) simply is 




(A2) 



i.e., we have a "cartesian" -like volume-element d£ dz for the normalization of $. 



With 



H{t) = K^ + K Z + V(t) 



(A3) 




(A4) 
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2 

m 



V(t) = - + — x + zE(t) sm(ut + if) (A6) 

the TDSE for z,t) assumes the form 

= (A7) 

The goal is to solve this TDSE. 

If A > 1/2 the transformation ( |A1| ) implies that $(0, z,t) — for all times. We discretize 
the (£, z)-space by 

& = jA£, j = l,2,...,J, z fe = (A; - K/2)Az, k = l,2,...,K (A8) 

with constant A£ and Az. While A = 1 yields the usual cylindrical coordinate system 
A = 3/2 turned out to offer the numerically more appropriate choice ||20|| . This is owing to 



the proper treatment of the wave function near the origin when the finite difference-formulas 
for the first and second derivatives in the Hamiltonian ( |A3| ) are applied to the wave function 
$. Note that uniform spacing in £ corresponds to non-uniform spacing in p. For A > 1 the 
p grid-width near the origin is smallest while it gets coarser far away from the origin. 

We use 3-point difference-formulas for all derivatives in K z and and impose as addi- 
tional boundary conditions 

= = HZ,ZK,t) = 0. (A9) 

In longer runs we apply a filter each time step which removes probability density moving 
towards the boundaries. This is a somewhat "shabby" method (similar to "imaginary po- 
tentials") but proper "absorbing boundary conditions" as discussed in [^] are not easily 
implemented in more than one dimensions. In any case, we always checked our numerical 
results upon sensitivity with respect to grid size and spacing. 

The time propagation is performed by applying the evolution operator 
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^ + Ai > - ^TAab^i G;SS9 11 " iAM( w)/2] (A10) 

with 

A(*)=Jf, + i^(t), 5(t) = ^ + ^(t) 

to the discretized representation of 2, t). This is the so-called Peaceman-Rachford 
method (PR) @, the alternating direction- version of the Crank-Nicholson-method for the 



TDSE in more than one dimension. The evolution operator (|A10|) is second order accurate 
in time and space (as long as the usual 3-point-difference formulas for the derivatives are 
used). Provided a non-iterative method for solving the implicit matrix equations 

(1 + iAtB{t n+1/2 )/2)^ n+l ' 2 = (1 - iAtA(t n+1/2 )/2)$ n , (All) 
(1 + iAL4(t„ +1/2 )/2)$ m+1 = (1 - iAtB(t n+1/2 )/2)^ 2 (A12) 

is chosen, the method is unconditionally stable. 

The stable ground state on our numerical grid was determined by propagating a "seed 
function" in imaginary time, i.e., we substituted At — > — iAt in ( |A10[ ). Here, renormalization 
of the wave function according (|A2|) after several time steps is necessary since imaginary time 
propagation is not unitary. Our experience was that during imaginary time propagation At 
had to be sufficiently small for $ converging to the ground state. A typical choice of our 
numerical parameters was (both for real and imaginary time propagation) 

A£ = Az = 0.1, At = 0.05, J = 60, K = 1000. 
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FIGURES 

FIG. 1. Ground state population T(t) vs time for an instantaneously switched on dc electric 
field. After a short transient behaviour (till ~ 2 atomic time units) the rates remain constant in 
time. The field strengths E as well as the constant rates W are indicated in the plot. 

FIG. 2. Contour plot of the probability density \ip(£, z)\ 2 after 15.5 atomic time units for 
the E = 0.3-case. The inlet shows the same situation as a surface plot. The electron escapes in 
negative z-direction by "over the barrier" -ionization. However, a peak remains at the Coulomb 
singularity. 

FIG. 3. The ground state populations in a strong ac field for two different peak field strengths 
{E = 0.3 and 0.5) and three different phases <p each. The dotted lines correspond to <p = 0, i.e., 
Esinuit, the dashed lines are the (p = 7r/2-case (Ecosut), and the intermediate case (p = 7r/4 is 
drawn dashed-dotted. 

FIG. 4. Contour plot of the probability density |VK£> z )\ 2 after 1 laser cycle for the 
E(t) = 0.3 cos wi-case. Owing to rescattered probability density wave packets have already formed. 
The inlet shows the corresponding surface plot of the probability density. 

FIG. 5. Ground state populations for hydrogen in a Gaussian laser pulse covering N cycles 
within the region where the electric field is 5% of the pulse amplitude E (see formulas (|l~7| ) and 
(HD for details). 

FIG. 6. Instantaneous ionozation rates vs the electric field present at the certain instant 
during the course of the laser pulse. The results have been obtained from different pulse shapes and 
frequencies: (+) sin 2 -pulse with oj = 0.2, (*) sin 2 -pulse with oj = 0.1, (O) instantaneously switched 
on dc field, (A) Gaussian pulse with to = 0.2. The curves are predictions from various analytical 
theories: (L) Landau, (Al) ADK, (A2) to BSI extended ADK, (K) Keldysh, (P) Posthumus, and 
(M) Mulser. The agreement in the region 0.15 < E < 0.5 is poor. The straight line is W = 2AE 2 
which fits the numerical data in this region quite well. 
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FIG. 7. Comparison of the numerically determined ground state populations vs time (drawn 
dotted) with the analytical predictions by means of the empirical formula ( f2l| ) (drawn solid). The 
dashed curve shows the result for the E = 0.5, N = 6-result when only the Landau rate (||) is 
applied during the entire pulse. 

FIG. 8. Comparison of the numerical square pulse-results with the predictions by formula 
(|2l|). In the upper plot (a) the agreement is very good while in the lower plot (b) formula (|2l| ) 
suffers from the transient ionization dynamics caused by the abrupt jump in the electric field at 
t = 0. 
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Fig. 1. D. Bauer and P. Mulser, "Exact field ionization rates ..." 
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Fig. 2. D. Bauer and P. Mulser, "Exact field ionization rates ..." 
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Fig. 3. D. Bauer and P. Mulser, "Exact field ionization rates ..." 
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Fig. 4. D. Bauer and P. Mulser, "Exact field ionization rates ..." 
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Fig. 6. D. Bauer and P. Mulser, "Exact field ionization rates ..." 
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Fig. 8. D. Bauer and P. Mulser, "Exact field ionization rates ..." 



